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ABSTRACT 

We study the central dark matter (DM) cusp evolution in cosmologically grown galactic halos. Nu- 
merical models with and without baryons (baryons+DM, hereafter BDM model, and pure DM, PDM 
model, respectively) are advanced from identical initial conditions, obtained using the Constrained 
Realization method. The DM cusp properties are contrasted by a direct comparison of pure DM 
and baryonic models. We find a divergent evolution between the PDM and BDM models within the 
inner fewxlO kpc region. The PDM model forms a i? -1 cusp as expected, while the DM in the BDM 
model forms a larger isothermal cusp R~ 2 instead. The isothermal cusp is stable until z ~ 1 when 
it gradually levels off. This leveling proceeds from inside out and the final density slope is shallower 
than —1 within the central 3 kpc (i.e., expected size of the cusp), tending to a flat core within 
~ 2 kpc. This effect cannot be explained by a finite resolution of our code which produces only a 5% 
difference between the gravitationally softened force and the exact Newtonian force of point masses 
at 1 kpc from the center. Neither is it related to the energy feedback from stellar evolution or angular 
momentum transfer from the bar. Instead it can be associated with the action of DM+baryon subha- 
los heating up the cusp region via dynamical friction and forcing the DM in the cusp to flow out and 
to 'cool' down. The process described here is not limited to low z and can be efficient at intermediate 
and even high z. 

Subject headings: cosmology: dark matter — galaxies: evolution — galaxies: formation — galaxies: 
halos — galaxies: interactions — galaxies: kinematics and dynamics 



1. INTRODUCTION 

Cosmological numerical simulations of pure dark mat- 
ter (DM) halo formation have led to "universal" den- 
sity profile which can be approximated by p(r) = 
Ap s /{r/R s ){l + r/R s ) 2 (Navarro et al. 1997, NFW), with 
R s being a characteristic "inner" radius where logarith- 
mic density slope is —2 and p s is the density at R s . This 
profile steepens outside i? s and tends to a cusp ~ 
toward the center, at i? C us P ~ 0.1i? s . While this pro- 
file remains invariant under mergers (El-Zant 2008), its 
origin is yet to be explained (e.g., Syer & White 1998; 
Nusser & Sheth 1999; Shapiro et al. 2004; Ascasibar et 
al. 2007) 

The NFW density profile is a source of an ongoing 
controversy — its universality in the numerical pure DM 
models is in apparent contradiction with at least some 
of the observations of disk galaxies and galaxy clusters, 
which exhibit rather flat density cores (e.g., Flores & Pri- 
mack 1994; Kravtsov et al. 1998; Salucci & Burkert 2000; 
Sand et al. 2002; de Blok et al. 2003; de Blok 2005, 2007; 
Simon et al. 2003; but see Rhee et al. 2004; answered by 
Gentile et al. 2004, 2005; de Naray et al. 2008). While 
a number of 'exotic' explanations involving less conven- 
tional physics have been proposed, attempts to resolve 
this discrepancy within the CDM cosmology have been 
made as well. The DM cusp leveling off was attributed 
to the DM-baryon interactions, such as a dynamical fric- 
tion of DM/baryon inhomogeneities (i.e., substructure) 
against the DM halo background (El-Zant et al. 2001, 
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2004; Tonini et al. 2006), stellar bar-DM interaction 
(Weinberg & Katz 2002; Holley-Bockelmann et al. 2005; 
but see Sellwood 2003; McMillan & Dehnen 2005; Du- 
binski et al. 2008), and baryon energy feedback (e.g., 
Mashchenko et al. 2006; Peirani et al. 2008). In this 
Letter, we test the first option (El-Zant et al.) — the 
effect of the baryon+DM substructure in the fully self- 
consistent numerical simulations of halo formation in the 
ACDM cosmology with WMAP3 paramaters (more de- 
tails in Romano-Diaz et al. 2008). 

Early arguments about change of the central density 
profile claimed that baryons drag DM in the so-called 
adiabatic contraction, steepening the DM density slope 
(e.g., Blumenthal et al. 1986). However, they have 
neglected the clumpy nature of the accreting material 
in the assembling halo, which can give rise to clump- 
background dynamical friction and energy transfer from 
the clumps to the background. The ability of baryons to 
radiate their internal energy increases the binding energy 
and acts as a 'glue' on the DM substructure (El-Zant et 
al. 2001). El-Zant et al. have used the Monte-Carlo tech- 
nique to calculate the leveling of the NFW cusp in the 
presence of such clumps and found this process to be ef- 
ficient over ~ 1 — 2 Gyr. The initial conditions presumed 
the existence of a NFW cusp, a spherically-symmetric 
DM halo and indestructible clumps. 

The NFW cusp is characterized by the "temperature" 
(i.e., DM dispersion velocities) inversion which makes it 
thermodynamically improbable but dynamically stable 
in the absence of energy transfer mechanism (El-Zant 
et al. 2001). Dynamical friction of sufficiently bound 
clumps can trigger such energy flux into the cusp, in- 
creasing its dispersion velocities and washing it out. El- 
Zant et al. found that the necessary condition for the 
dynamical friction to have an effect is the aggregation of 
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Fig. 1.— Redshift evolution of DM density profiles, p(R) (left), p(R)R (middle) and p(R)R? (right) in PDM and BDM models: z = 3.55 
(solid), 2.12 (dotted), 1.0 (dashed), 0.61 (dot-dashed), 025 (dot-dash-dotted) and (long dashed). The PDM and BDM curves are displaced 
vertically for clarity. The inner 40 kpc of halos are shown. The vertical coordinate units are logarithmic and arbitrary. For the PDM 
model, the density is well fitted by the NFW profile over a large range in z, and R s ~ 28 kpc at z = 0. For the BDM model, the NFW fit 
is worse and i?i so ~ 15 kpc at the end. The insert provides p within 200 kpc range for a comparison. 



mass into the clumps more massive than ~ 10 of the 
DM halo — result confirmed for the galaxy clusters as 
well (El-Zant et al. 2004). 

In fully self-consistent numerical simulations with 
baryons, the question is whether the NFW cusp forms in 
the first place. Such simulations have been either limited 
to high redshifts (z — 3.3, Gnedin et al. 2004) or focused 
on the adiabatic contraction and had an insufficient res- 
olution (Gustafsson et al. 2006). They have verified that 
baryons lead to a steeper DM density profile than the 
NFW cusp, but left open its subsequent evolution. This 
issue is addressed here. 

2. NUMERICS AND INITIAL CONDITIONS 

Numerical simulations have been performed with the 
FTM 4.5 hybrid iV-body/SPH code (Heller & Shlosman 
1994; Heller et al. 2007; Romano-Diaz et al. 2008) us- 
ing physical and not comoving coordinates. The num- 
ber of DM particles is 2.2 x 10 6 and the SPH particles 
- 4 x 10 5 . The gravity is computed with the falcON 
routine (Dehnen 2002) which scales as O(N). The grav- 
itational softening (PI of falcON) is e gra v = 500 pc, for 
DM, stars and gas. The calculated force differs by ~ 5% 
from the exact Newtonian force between point masses at 
R ~ 2e grav = 1 kpc. We assume the ACDM cosmology 
with WMAP3 parameters, Q m = 0.24, Q A = 0.76 and 
h = 0.73, where h is the Hubble constant in units of 
100 km s -1 Mpc _1 . The variance cr§ = 0.76 of the den- 
sity field convolved with the top hat window of radius 
8ft." 1 Mpc -1 was used to normalize the power spectrum. 
The star formation algorithm is described in Heller et 
al. (2007). Several generations of stars are allowed to 
form from an SPH particle. The energy and momentum 
feedback into the ISM is implemented, and the relevant 
parameters are (Heller et al.): the energy thermalization 
6sf = 0.3, the cloud gravitational collapse ag = 1, and 
the self-gravity fudge factor a cr i t = 0.5. 

The initial conditions generated here are those of 
Romano-Diaz et al. (2008): we use the Constrained Re- 
alizations (CR) method (Hoffman & Ribak 1991) within 
a restricted volume of 8 3 /i _1 Mpc, where a 5/i _1 Mpc 
sphere is carved out and evolved. The constructed Gaus- 
sian field is required to obey a set of constraints of ar- 
bitrary amplitudes and positions (e.g., Romano-Diaz et 



al. 2006, 2007). Two constraints were imposed on the 
initial density field, first — that the linear field Gaus- 
sian smoothed with a kernel of 1.0 x 10 12 /i _1 M a has an 
over-density of S = 3 at the origin (2.5cr perturbation, 
where a 2 is the variance of the appropriately smoothed 
field). It was imposed on a 256 3 grid and predicted to 
collapse by z c ~ 1.33 by the top- hat model. This pertur- 
bation is embedded in a region (2nd constraint) corre- 
sponding to 5 x 10 13 ft _1 M in which the over-density is 
zero, i.e., the unperturbed universe. The random com- 
ponent of the CR introduces density perturbations on 
all scales, thus leading to major mergers (Romano-Diaz 
et al. 2006, 2007). The mass inside the computational 
sphere is ~ 6.1 x 10 12 h^M®. In the baryonic model, 
we have randomly replaced 1/6 of DM particles by equal 
mass SPH particles. Hence, f2 m is not affected. 

3. DM DENSITY AND VELOCITY DISPERSION PROFILES 

Two models analyzed here are those of a pure DM 
(hereafter PDM) and DM+baryons (BDM). Evolution 
of DM density profiles for both models is shown in Fig. 1 
at various redshifts and in three different ways: as p(R) 
per se, factored by R, and by R 2 , in order to emphasize 
the and R~ 2 cusps. At early redshifts, z ~ 7, the 
DM density profiles of both models are nearly identical. 
With time, the BDM halo becomes more cuspy — the 
DM experiencing the adiabatic contraction by baryons. 
This region of higher density extends gradually to larger 
R, up to ~ 15 kpc at z — 0. Between z ~ 4 and 1 the 
DM density profile becomes and remains nearly isother- 
mal within this region — p(R)R 2 is flat there (Fig. 1, 
right frame). In comparison, the PDM halo exhibits the 
slope of — 2 at i? s ~ 28 kpc only, as expected, with the 
NFW cusp size R cusp ~ 3 kpc at z = 0. Hence over 
prolonged time period, while the PDM model displays 
the NFW cusp, the DM in the BDM model remarkably 
forms an isothermal cusp, down to the central kpc, i.e., 
to ~ 2e grav . The isothermal cusp in the BDM model is 
gradually erased inside i?j so ~ 15 kpc after z ~ 1. In fact, 
within the central 3 kpc it becomes flatter than — 1, i.e., 
flatter than the NFW cusp, as seen in all three frames 
of Fig. 1, and within ~ 2 kpc forms a rather flat core. 
At z = 0, the PDM model is denser outside i?i so , while 
BDM remains denser within this radius. 
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Fig. 3. — Redshift evolution of the power index |/9| in the PDM 
(blue solid) and BDM (dashed red) models. /3 is calculated by 
fitting a power law o"p M ~ R~@ between 1 kpc and 10 kpc. Note 
that i§ > for BDM and /3 < for PDM models. 



Fig. 2. — Redshift evolution of DM velocity dispersions in PDM 
(left) and BDM (right) models. Except for the lowest ones, the 
curves are displaced vertically up for clarity. The second curves 
(from the bottom) are displaced by a factor of 2, the third — by 
a factor of 2 2 , the fourth — by a factor of 2 3 , and the last ones 
— by a factor of 2 4 . The colored width represents a la dispersion 
around the mean. The inner 200 kpc of halos are shown. The 
vertical coordinate units are logarithmic. 



Furthermore, in the PDM model, the density around 
1 kpc from the center stays virtually unchanged after 
z ~ 3. The baryonic model shows a somewhat different 
behavior. Namely, it keeps the same DM density at 1 kpc 
between z ~ 3 — 0.7, which decays thereafter by a factor 
of - 3. 

Behavior of the velocity dispersions, <tdm, in both 
models mirrors that of the density. The NFW DM cusp 
in the PDM model forms early and is characterized by 
the 'temperature' inversion as shown in Fig. 2, where 
we display the z evolution of the DM dispersion veloci- 
ties (see also El-Zant et al. 2001). In contrast, the BDM 
model shows a cuspy distribution of <tt>m(R), which after 
z — 4 can be approximated well by a single power law, 
(Jum(R) 2 ~ R 13 inside central 10 kpc. The evolution 
of |/3 j is shown in Fig. 3 for both models. It is steadily 
increasing in the PDM model until the end of the major 
mergers epoch at z — 1.5, where it levels off. The BDM 
model shows a similar increase in |/3|, but well beyond 
the mergers epoch, until z — 0.7, where it sharply de- 
clines. This decline in /3 is crucial in understanding the 
fate of the density cusp, as we discuss below. 

4. DISCUSSION 

We have performed a direct comparison between the 
DM halo evolution in the pure DM and baryonic models, 
from identical cosmological initial conditions. Our goal is 
to understand the effect of the baryons on the DM den- 
sity profile. We, therefore, focus on the DM evolution 
within the central region, quantify the cuspiness of the 
density distribution there, and attempt to relate its evo- 
lution to that of the DM substructure. Our results show 
that the DM in the PDM model has established a NFW 
cusp early, and the only subsequent changes we observe 
are related to sudden increases in R s at the time of the 



major mergers (e.g., Romano-Diaz et al. 2006, 2007). On 
the other hand, the baryonic (BDM) model goes through 
a two-step process: first, in the central — 10 kpc, the DM 
is dragged in by the baryons in what fits the adiabatic 
contraction (e.g., Gnedin et al. 2004; Gustafsson et al. 
2006). This stage is concurrent with the formation and 
growth of a galactic disk, but we did not verify whether 
the disk is solely responsible for the DM contraction. For 
an extended period of time, after z — 4, an isothermal 
cusp forms with pdm ~ R 2 and persists. In the next 
step, the DM cusp is being gradually washed out from 
inside out, mainly after z — 1. By z = 0, the isothermal 
cusp is largely erased, and the density slope within the 
central 10 kpc is less steep than —2. In fact within the 
central — 3 kpc, it is less steep than the NFW slope of —1. 
We note, that the DM in the BDM model does not form 
the NFW cusp in the first place, but rather by-passes it 
in favor of an isothermal DM cusp which is subsequently 
washed out. For example, the 'temperature' inversion 
characteristic of the NFW cusp is never observed in the 
BDM model (Fig. 2). 

We first argue, that this leveling of the DM cusp is 
not a numerical artifact, i.e., is not caused by a finite 
resolution of the code which is 2e grav = 1 kpc. At this 
distance, softened and exact Newtonian forces between 
point masses differ by 5% — this cannot account for the 
observed changes in density profile further out. We also 
note that the DM density peak in the BDM model is 
occasionally offset from the baryonic peak by a few kpc, 
around z *~ 1 and especially 0.4, with the subhalos in- 
flux into the center, in what can be treated as m = 1 
instability. Therefore, we have fitted the DM halo pro- 
file based on the position of the DM density peak at 
each frame. The cusp will be 'smeared' if this procedure 
is not followed (McMillan & Dehnen 2005). Also, the 
prime halo moves with respect to the center of mass of 
the computational sphere, this is related to the residual 
large-scale streaming motions in the DM (Romano-Diaz 
et al. 2008). 

Next, we note that most of the cusp leveling happens 
after z — 1 (Fig. 1). This correlates with a decrease in 
the slope of the central dispersion velocities (Fig. 2). The 
slope change in [3 (Fig. 3) comes from the decrease in the 
central dispersion velocities around few kpc, while ctdm 
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Fig. 4. — Evolution of the subhalo number iV s bh within central 
30 kpc of the prime halo with redshift — PDM (lower blue) and 
BDM (upper red). Note the sharp increase in AT s bh for the BDM 
model after z ~ 1, 0.5 and 0.2. 



at R ~ 10 kpc remains unchanged. What is the reason 
for this apparent 'cooling' of the isothermal DM cusp? 

Any decrease in central ctdm of a dissipationless self- 
gravitationg entity in a virial equilibrium is counter- 
intuitive. It clearly has nothing to do with the major 
mergers, as this epoch ends at z ~ 1.5 in our models. 
Any subsequent smooth accretion of dissipative baryons 
will only increase cdm as a result of the ongoing adia- 
batic contraction. However, if one considers accretion of 
clumpy baryons, which act as a 'glue' on the DM subha- 
los where they reside, additional process of 'heating' the 
DM by means of a dynamical friction becomes important. 
The heating rate is skewed to the density peak and will 
cause the DM to stream out of the cusp's gravitational 
well, eventually decreasing the 'temperature' gradient at 
the center. Fig. 3 displays this behavior of j3 after it 
reaches its peak at z ~ 0.7. For a comparison, the PDM 
halo shows a flat j3 after the epoch of major mergers, for 
nearly 10 Gyr. 

Central region crossing by subhalos (as well as accre- 
tion and tidal disruption) continues to the present time. 
Why do changes in the cusp become significant after 
2 ~ 1? Fig. 4 displays the evolution of the subhalo pop- 
ulation within the central 30 kpc of the BDM halo. We 
observe that at least in two instances, z ~ 1 and 0.4, 
there is a clear excess in the BDM subhalo population, 
corresponding to the splash in the subhalo influx rates. 
The number of these subhalos also exceeds that of the 
PDM by ~ 2. The BDM subhalos appear in waves within 
the central region and these waves coinside with the DM 
streaming out of the center and 'cooling' down. Hence, 
the process of washing out the DM cusp correlates with 
the influx of subhalos into the innermost region of the 
BDM halo. We show elsewhere that waves of subhalos 
crossing the central region of the prime halo originate 
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in the filament, cluster there and enter the prime halo 
before they merge among themselves. 

We find that the mass distribution of subhalos (in the 
computational sphere) with masses M s bh evolves with 
redshift — at high z J> 4 it can be approximated by a 
power law, N shh ~ Af"*, both in PDM and BDM models 
(Romano-Diaz et al. 2008). But in the central virialized 
100 kpc of the halo it evolves differently from the field 
subhalos. Specifically, it is heavily skewed toward more 
massive ones after z ~ 1, while that of the PDM subha- 
los is only weakly so. This is significant, because the effi- 
ciency of dynamical friction increases substantially with 
the clump-to-primc halo virial mass ratio (El-Zant et al. 
2001). At z ~ 4, this ratio for the most massive subhalo 
residing within the virialized halo is ~ 0.06. By z ~ 1 
this ratio is ~ 10~ 3 , and at z = it is ~ 8 x 10 — still 
higher than the minimum required by El-Zant et al. 

The size of the obtained core in the DM distribution of 
the BDM model can be compared with observationally 
inferred cores. Following Salucci et al. (2007), for the 
virial mass of our halo at z = 0, M v i r ~ 4 x 10 12 Mq, 
the expected core is ~ 1.6 kpc. Spano et al. (2008) lists 
a large range of core sizes, starting from ~ 1 kpc. Both 
are compatible with the value obtained here. 

In summary, we find that the DM density cusp corre- 
sponding to -Rcusp in the PDM model is leveled off by the 
action of subhalos in the presence of baryons in the BDM 
model. The energy feedback from stellar evolution has 
decayed by this time without any effect on the cusp. We 
do not detect the angular momentum (J) transfer from 
the bar to the cusp — the latter J stays constant in 
time. We also find that the BDM model forms a steeper 
isothermal rather than NFW-type cusp for extended time 
period. It goes through a two step process which involves 
gravitational contraction, followed by the influx of sub- 
halos which heat up and dissolve the cusp. 

Our models, although containing about 1.6 x 10 6 par- 
ticles per halo, are still insufficient to fully resolve the 
formation of substructure around the prime halo. They 
should be viewed rather as a lower limit on the efficiency 
of the process reported here. Neither is the process de- 
scribed here limited to low z one can easily envision 
scenario (s) when the subhalos act at intermediate and 
even high z. Overall, this work stands in agreement with 
Gnedin et al. (2004) who stopped the simulation at red- 
shift z = 3.3, well before the leveling of the density cusp. 
However, Gustafsson et al. (2006) barely lacked reso- 
lution to detect the cusp behavior observed here. We 
have recently learned that J. P. Ostriker (priv. comm.) 
obtained similar results on the cusp flattening. 
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